      PROGRAM AUTOFLAT
C Vyvoreni autoflat-fieldu z mnoha CCD snimku.
C
C FH, Rijen 1997
C FH, duben 1998
C FH, cerven 1998
C
C Porad je to jen experimenatlni verze!
! Opatrnost nutna!
!
     
       IMPLICIT NONE

       INTEGER MAXX,MAXY,NSIDE
       PARAMETER( MAXX = 2000, MAXY = 2000, NSIDE = 100 )

       CHARACTER*80 SB,BUF
       LOGICAL SIMPLE,EXTEND,ANYF,EXIST,FLAGVAL
       INTEGER I,J,L,K,N,LN,II,JJ
       INTEGER STATUS,BITPIX,NAXIS,NAXES(2) ! dimenze snimku
       INTEGER IMAGE(MAXX,MAXY),IBUF(MAXX*MAXY)
       INTEGER BLOCKSIZE,PCOUNT,GCOUNT,FPIXELS(2),LPIXELS(2)
       INTEGER NAXISX,NAXESX(2),VALUE,MED,MEDIAN
       REAL BSCALE,BZERO,FFLAT(MAXX,MAXY)
       DOUBLE PRECISION S,DS,DS1
       
       EXTERNAL MEDIAN

c       DATA FFLAT/4 000 000*0.0/
       DATA NAXES/0,0/

       DO J = 1,MAXY
         DO I = 1,MAXX
           FFLAT(I,J) = 0.0
         ENDDO
       ENDDO

C otevreni souboru se jmeny souboru na vyrobu flat-fieldu
c       OPEN(1,FILE='af.in',STATUS='OLD')

       N = 0
! nacteni jmena souboru
   10  READ(*,'(A)',END=20) SB   
         if( sb .eq. ' ' ) goto 10 
         print *,sb(:50)
! otevreni souboru
         CALL FTOPEN(2,SB,0,BLOCKSIZE,STATUS)
         IF( STATUS.NE.0 )THEN
           PRINT *,'STATUS41=',STATUS
           STOP
         ENDIF    
! cteni parametru
         CALL FTGHPR(2,2,
     &        SIMPLE,BITPIX,NAXISX,NAXESX,PCOUNT,GCOUNT,EXTEND,STATUS)
         IF( STATUS.NE.0 )THEN
            PRINT *,'Spatne dimenze obrazku.'
            PRINT *,'STATUS48=',STATUS
            STOP
         ENDIF
         CALL FTGKYE(2,'BSCALE',BSCALE,BUF,STATUS)
         CALL FTGKYE(2,'BZERO',BZERO,BUF,STATUS)
! cteni obrazku
         L = 0
         DO J = 1,NAXESX(2)
           DO I = 1,NAXESX(1)
             L = L + 1
             CALL FTGPVJ(2,1,L,1,0,VALUE,ANYF,STATUS)
             IF( STATUS.NE. 0 .OR. ANYF )THEN
               PRINT *,'STATUS=',STATUS
c               STOP
             ENDIF             
             IMAGE(I,J) = VALUE
           ENDDO
         ENDDO

      DO J = 1,NAXESX(2)
        DO I = 1,NAXESX(1)
          IF( I/(NSIDE/2)+1.NE.I .OR. J/(NSIDE/2)+1.NE.J ) GOTO 60   
! pro kazdy 100-ty pixel obrazku se udela lokalni median
          L = 0
          DO JJ = -NSIDE,NSIDE
            DO II = -NSIDE,NSIDE
              IF(   J+JJ.GE.1 .AND. J+JJ.LE.NAXESX(2) )THEN
                IF( I+II.GE.1 .AND. I+II.LE.NAXESX(1)   )THEN
                  L = L + 1
                  IBUF(L) = IMAGE(J+JJ,I+II)
                ENDIF
              ENDIF
            ENDDO
         ENDDO 
         MED = MEDIAN(L,IBUF,L/2)
! iteracne spocte lokalni odchylku
         DS1 = 1e10
         ds = ds1
 50      S = 0.0
         LN = 0
         DO JJ = -NSIDE,NSIDE
           DO II = -NSIDE,NSIDE
              IF(   J+JJ.GE.1 .AND. J+JJ.LE.NAXESX(2) )THEN
                IF( I+II.GE.1 .AND. I+II.LE.NAXESX(1)   )THEN
                  L  = ABS(IMAGE(I+II,J+JJ) - MED)
                  IF( L .LT.5.0*DS ) THEN
                    S = S + L
                    LN = LN + 1
                  ENDIF
                ENDIF
              ENDIF
           ENDDO
         ENDDO   
         DS = S/LN
         s = abs(ds - ds1)
         ds1 = ds
         IF( s/DS.GT.0.2 ) GOTO 50
! a hodnotu v obrazku pricte k formovanemu flat-fieldu
 60      VALUE = IMAGE(I,J)
         IF( ABS(VALUE - MED).GT.5.0*DS ) VALUE = MED
         FFLAT(I,J) = FFLAT(I,J) + REAL(VALUE)/REAL(MED)
        ENDDO
      ENDDO
      N = N + 1
      CALL FTCLOS(2,STATUS)
      GOTO 10
   20 CONTINUE 
c       print *,i,j,fflat(i,j)
       DO J = 1,NAXESX(2)
         DO I = 1,NAXESX(1)
            IMAGE(I,J) = REAL(MED*FFLAT(I,J))/REAL(N)
         ENDDO
       ENDDO

       INQUIRE(FILE='autoflat.fts',EXIST=EXIST)
       IF( EXIST ) THEN
         OPEN(9,FILE='autoflat.fts')
         CLOSE(9,STATUS='DELETE')
       ENDIF
C Vytvoreni vystupniho obrazku
       CALL FTINIT(9,'autoflat.fts',1,STATUS)
       IF( STATUS.NE. 0 )THEN
         PRINT *,'STATUS103=',STATUS
         STOP
       ENDIF
C zapsani informace do hlavicky
       SIMPLE = .TRUE.
       BITPIX = 16
       NAXIS = 2
       NAXES(1) = NAXESX(1)
       NAXES(2) = NAXESX(2)
       CALL FTPHPS(9,BITPIX,NAXIS,NAXES,STATUS)
       CALL FTPKYS(9,'OBJECT','FLAT-FIELD','',STATUS)
       CALL FTPKYE(9,'BSCALE',BSCALE,10,'',STATUS)
       CALL FTPKYE(9,'BZERO',BZERO,10,'',STATUS)

C zapsani  dat
       CALL FTP2DJ(9,1,MAXX,NAXES(1),NAXES(2),IMAGE,STATUS)
       IF( STATUS.NE. 0 )THEN
         PRINT *,'STATUS120=',STATUS
         STOP
       ENDIF
C Zavreni souboru
       CALL FTCLOS(9,STATUS)

       print *,'Bye.'
       END       


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      INTEGER FUNCTION MEDIAN(N,A,K)
C find Kth element from array A
C According to procedure in
C N.Wirth: Algorithms + Data Structures = Programs, Prentice-Hall, 1975
C
      INTEGER N,K,L,R,I,J,W,X,A(*)

      L = 1
      R = N
10    IF( L.LT.R )THEN
        X = A(K)
        I = L
        J = R
30      IF( A(I).LT.X ) THEN
            I = I+1
            GOTO 30
          ENDIF
40        IF( X.LT.A(J) )THEN
            J = J-1
            GOTO 40
          ENDIF
          IF( I.LE.J )THEN
            W = A(I)
            A(I) = A(J)
            A(J) = W
            I = I+1
            J = J-1
          ENDIF
        IF( I.LE.J ) GOTO 30
        IF( J.LT.K ) L = I
        IF( K.LT.I ) R = J
        GOTO 10
      ENDIF
      MEDIAN = X
      END
C

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC



